Data Clinic 4: Image processing

PART A: Image basics, histograms

1. Reading the image and displaying it

We will start off by reading our images and displaying them in full color.

library(EBImage)
library(CRImage)

setwd("/Volumes/SD/GoogleDrive/0.DA/Assignments/4")
image <- readImage('image1.jpg')
myImage <- readImage('cat.jpg')

display(image, method="raster")

display(myImage, method="raster")

2. Converting the colored image to gray values and display it

Now we will convert the images to grayscale and display them.

grayImage <- EBImage::channel(image,"gray")
myGrayImage <- EBImage::channel(myImage,"gray")
#If there is masking of channel from other packages use the following one:
#grayImage <- EBImage::channel(image,"gray")
display(grayImage, method="raster")

display(myGrayImage, method="raster")

Compared to our colored image which had dimensions of 640 x 427 x 3, our new image has dimensions 640 x 427 x 1.
This is because for each pixel we had 3 values (RGB), but now we have only 1 value (grayscale).

3. Converting the grayscale image to different binary images

We will convert our grayscale images to binary images using thresholds.
For this we will first create a function that does this.

#convert to Binary
threshold_values <- function(image, threshold){
  image[image>threshold]=1
  image[image<=threshold]=0
  return (image)
}
bin02=threshold_values(grayImage,0.2)

bin04=threshold_values(grayImage,0.4)

bin06=threshold_values(grayImage,0.6)

bin08=threshold_values(grayImage,0.8)

bin09=threshold_values(grayImage,0.9)

par(mfrow = c(5,2))
plot(bin02)
title("Threshold = 0.2")
plot(image)
title(main = "Original")

plot(bin04)
title("Threshold = 0.4")
plot(image)
title(main = "Original")

plot(bin06)
title("Threshold = 0.6")
plot(image)
title(main = "Original")

plot(bin08)
title("Threshold = 0.7")
plot(image)
title(main = "Original")

plot(bin09)
title("Threshold = 0.9")
plot(image)
title(main = "Original")

bin02=threshold_values(myGrayImage,0.2)

bin04=threshold_values(myGrayImage,0.4)

bin06=threshold_values(myGrayImage,0.6)

bin08=threshold_values(myGrayImage,0.8)

bin09=threshold_values(myGrayImage,0.9)

par(mfrow = c(5,2))
plot(bin02)
title("Threshold = 0.2")
plot(myImage)
title(main = "Original")

plot(bin04)
title("Threshold = 0.4")
plot(myImage)
title(main = "Original")

plot(bin06)
title("Threshold = 0.6")
plot(myImage)
title(main = "Original")

plot(bin08)
title("Threshold = 0.7")
plot(myImage)
title(main = "Original")

plot(bin09)
title("Threshold = 0.9")
plot(myImage)
title(main = "Original")

4. Contrast Stretching, Clipping and Range Compression

Now we will apply some digital image processing functions.

Contrast Stretching

We will first define a function that does contrast stretching:

contrast_stretching <- function(image,a,b,y_a,y_b,L, alpha, beta, gamma){
  a = a / L
  b = b / L
  y_a = y_a / L
  y_b = y_b / L
  if(a > 1 || b > 1){
    print("Error")
  }
  image[image < a] <- alpha * image[image < a]
  image[image >= a & image < b ] <- beta * (image[image >= a & image < b ] - a) + y_a 
  image[image >= b] <- gamma * (image[image >= b] - b) + y_b
  return(image)
}

Now we will apply this to our images:

contrast_stretched_img1 <- contrast_stretching(grayImage, a = 50, b = 150, y_a = 30, y_b = 200,L = 255, alpha=0.5, beta=5, gamma=20) 
par(mfrow = c(2,2))
display(contrast_stretched_img1, method = "raster")
title("Contrasted")
display(grayImage, method = "raster")
title(main = "Original")


contrast_stretched_img1 <- contrast_stretching(myGrayImage, a = 50, b = 150, y_a = 30, y_b = 200,L = 255, alpha=0.5, beta=5, gamma=20) 
par(mfrow = c(2,2))

display(contrast_stretched_img1, method = "raster")
title("Contrasted")
display(myGrayImage, method = "raster")
title(main = "Original")

Clipping

We are first going to define the function that performs clipping:

clipping <- function(image,a,b,L,beta){
  a = a / L
  b = b / L
  if(a > 1 || b > 1){
    print("Error")
  }
  image[image < a] <- 0
  image[image >= a & image < b ] <- beta * (image[image >= a & image < b ] - a)
  image[image >= b] <- beta * (b - a )
  return(image)
}

Now we can apply it to our images:

clipped_img1 <- clipping(grayImage, a = 0.4, b = 0.6, L = 1.0, beta=2) 
par(mfrow = c(2,2))
display(clipped_img1, method = "raster")
title("Clipped")
display(grayImage, method = "raster")
title(main = "Original")

clipped_img1 <- clipping(myGrayImage, a = 0.4, b = 0.6, L = 1.0, beta=2) 
par(mfrow = c(2,2))

display(clipped_img1, method = "raster")
title("Clipped")
display(myGrayImage, method = "raster")
title(main = "Original")

Compression

We first define the function that does compression:

compress <- function(image,c,L){
  image <- c * log10(1 + image)
  return(image)
}

Now we can apply it to our images:

compressed_img1 <- compress(grayImage, c = 50, L = 1.0) 
par(mfrow = c(2,2))
display(compressed_img1, method = "raster")
title("Compressed")
display(grayImage, method = "raster")
title(main = "Original")

compressed_img1 <- compress(myGrayImage, c = 5, L = 1.0) 
par(mfrow = c(2,2))

display(compressed_img1, method = "raster")
title("Compressed")
display(myGrayImage, method = "raster")
title(main = "Original")

We used a much lower value of c for our image because it was overall brighter. Very high values of c made the picture completely white otherwise.

5. Histogram equalization

We will now perform histogram equalization on our images.

First let’s plot the normal and equalized histogram for the original picture.

print('Image - histogram')

## [1] "Image - histogram"

equalized_grayImage = equalize(grayImage, range = c(0, 1), levels = 256)
par(mfrow = c(2,2))
hist(grayImage)

hist(equalized_grayImage)

plot(grayImage)
title("original")

plot(equalized_grayImage)
title("equalized")

Now we are going to see the same histograme for our image.

print('Our image - histogram')

## [1] "Our image - histogram"

equalized_grayImage = equalize(myGrayImage, range = c(0, 1), levels = 256)
par(mfrow = c(2,2))
hist(myGrayImage)
hist(equalized_grayImage)

plot(myGrayImage)
title("original")

plot(equalized_grayImage)
title("equalized")

We can notice that histogram equivalization lowered the brightness on our picture.

PART B: Filters, Convolution

We will use a new picture for the this part.

image <- readImage('image2.jpg')
grayImage <- EBImage::channel(image,"gray")
display(x = image, method="raster")

6. Filtering

We will visualize intensities across positions in a row of our matrix.
Then we will apply a filter to this and see how it looks.

First let’s start with the original picture:

row <- grayImage[200,]
par(mfrow = c(2,1))
plot(seq(length(row)),row, xlab='Position', ylab='Intensity', type="l")
title("Original")
filter9 <- matrix( c(1/9) ,nrow=1,ncol = 9)
row_filtered_9 <- filter2(x = row, filter = filter9)
plot(seq(length(row)),row_filtered_9, xlab='Position', ylab='Intensity', type="l")
title("Filtered")

par(mfrow = c(2,1))
row <- myGrayImage[200,]
plot(seq(length(row)),row, xlab='Position', ylab='Intensity', type="l")
title("Original")

filter9 <- matrix( c(1/9) ,nrow=1,ncol = 9)
row_filtered_9 <- filter2(x = row, filter = filter9)
plot(seq(length(row)),row_filtered_9, xlab='Position', ylab='Intensity', type="l")
title("Filtered")

We can clearly observe how filtering reduces the small variations and smoothes out the line.

7. Blur Filtering

We will apply a 2D blur filter that will average out the neighouring pixels.

par(mfrow = c(1,2))
filter2d <- matrix( c(1/81) ,nrow=9,ncol = 9)
grayImage_2d <- filter2(x = grayImage, filter = filter2d)
display(grayImage_2d, method = "raster")
title('Original image blurred')

filter2d <- matrix( c(1/81) ,nrow=9,ncol = 9)
grayImage_2d <- filter2(x = myGrayImage, filter = filter2d)
display(grayImage_2d, method = "raster")
title('Our image blurred')

Now we will apply a 1D Filter whose values sum to 1, and apply again the same filter but transposed to the convolved image from the first filter.

par(mfrow = c(1,2))
filter1d <- matrix( c(1/9), nrow=9, ncol = 1)
grayImage_1d10 <- filter2(x = grayImage, filter = filter1d)
display(grayImage_1d10, method = "raster")
title("1D filter on original image")
grayImage_1d10_transp <- filter2(x = grayImage_1d10, filter = t(filter1d))
display(grayImage_1d10_transp, method = "raster")
title("Extra 1D filter on previous (convolved) image")

We can notice that after the two filters the image is very similar to the image produced by the 2D filter.
But after doing some precise checks in R it doesn’t seem like they are the same.

par(mfrow = c(1,2))
filter1d <- matrix( c(1/9), nrow=1, ncol = 9)
grayImage_1d10 <- filter2(x = myGrayImage, filter = filter1d)
display(grayImage_1d10, method = "raster")
title("1D filter on original image")
grayImage_1d10_transp <- filter2(x = grayImage_1d10, filter = t(filter1d))
display(grayImage_1d10_transp, method = "raster")
title("Additional 1D filter on previous (convolved) image")

8. High pass filtering

We will now apply a high pass filter to our data.
But after applying it, we will normalize the data back into a normal range (0 to 1).

par(mfrow = c(1,2))
hp_filter <- matrix(c(0,-1,0,-1,4,-1,0,-1,0),nrow=3,ncol=3)
hp_image <- filter2(grayImage, hp_filter)
min <- min(hp_image)
max <- max(hp_image)
hp_image <- (hp_image - min) / (max - min)

display(grayImage, method="raster")
title("Original image")

display(hp_image, method="raster")
title("High pass filtered image")

par(mfrow = c(1,2))
hp_filter <- matrix(c(0,-1,0,-1,4,-1,0,-1,0),nrow=3,ncol=3)
hp_image <- filter2(myGrayImage, hp_filter)
min <- min(hp_image)
max <- max(hp_image)
hp_image <- (hp_image - min) / (max - min)

display(myGrayImage, method="raster")
title("Original image")

display(hp_image, method="raster")
title("High pass filtered image")

The filter we applied is the presents the edges in the image or more precisely places where the (second derivative of) change of intensity is large.

PART C: Edge detection

We will use a new picture for the this part.

image <- readImage('image3.png')
grayImage <- EBImage::channel(image,"gray")
display(x = image, method="raster")

9. Calculating gradient magnitude and showing it for various thresholds

We will calculate the gradient magnitude after computing the local gradients.
To compute the local gradients we will use simple 1D arrays of (-1,0,1) and (1,0,-1) as shown on lecture slides.

filter_x <- matrix(c(-1,0,1),nrow=1,ncol=3)
filter_y <- matrix(c(1,0,-1),nrow=3,ncol=1)

gradient_x <- filter2(x = grayImage, filter = filter_x)
gradient_y <- filter2(x = grayImage, filter = filter_y)

gradient_mag <- sqrt(gradient_x^2+gradient_y^2)

#display(x = gradient_x, method="raster")
#display(x = gradient_y, method="raster")
print("Gradient magnitudes of Lena")

## [1] "Gradient magnitudes of Lena"

display(x = gradient_mag ,method="raster")
title("Lena")

Now let’s explore the gradient magnitude at various thresholds:

First we start with Lena:

grad_mag_1 <- threshold_values(gradient_mag,0.05)

grad_mag_2 <- threshold_values(gradient_mag,0.25)

grad_mag_3 <- threshold_values(gradient_mag,0.5)

grad_mag_4 <- threshold_values(gradient_mag,0.75)

grad_mag_5 <- threshold_values(gradient_mag,1.0)

par(mfrow = c(5,2))
plot(grad_mag_1)
title("Threshold = 0.05")
plot(grayImage)
title(main = "Original")

plot(grad_mag_2)
title("Threshold = 0.25")
plot(grayImage)
title(main = "Original")

plot(grad_mag_3)
title("Threshold = 0.5")
plot(grayImage)
title(main = "Original")

plot(grad_mag_4)
title("Threshold = 0.75")
plot(grayImage)
title(main = "Original")

plot(grad_mag_5)
title("Threshold = 1.0")
plot(grayImage)
title(main = "Original")

We can notice that only stronger (intensive) edges show up when we apply higher thresholds.

Now we will do the same for our image:

filter_x <- matrix(c(-1,0,1),nrow=1,ncol=3)
filter_y <- matrix(c(1,0,-1),nrow=3,ncol=1)

gradient_x <- filter2(x = myGrayImage, filter = filter_x)
gradient_y <- filter2(x = myGrayImage, filter = filter_y)

gradient_mag <- sqrt(gradient_x^2+gradient_y^2)

display(x = gradient_mag ,method="raster")
title("My image")

grad_mag_1 <- threshold_values(gradient_mag,0.05)

grad_mag_2 <- threshold_values(gradient_mag,0.25)

grad_mag_3 <- threshold_values(gradient_mag,0.5)

grad_mag_4 <- threshold_values(gradient_mag,0.75)

grad_mag_5 <- threshold_values(gradient_mag,1.0)

par(mfrow = c(5,2))
plot(grad_mag_1)
title("Threshold = 0.05")
plot(myGrayImage)
title(main = "Original")

plot(grad_mag_2)
title("Threshold = 0.25")
plot(myGrayImage)
title(main = "Original")

plot(grad_mag_3)
title("Threshold = 0.5")
plot(myGrayImage)
title(main = "Original")

plot(grad_mag_4)
title("Threshold = 0.75")
plot(myGrayImage)
title(main = "Original")

plot(grad_mag_5)
title("Threshold = 1.0")
plot(myGrayImage)
title(main = "Original")

The value of 0.25 as a threshold seems rather good for both images.
Higher thresholds are not very informative because gradient maginuted rarely have values above those thresholds.

2 Function with Prewitt and Sobel operators

Now we will implement a function that allows the user to choose between Prewitt and Sobel operators, and allows the user to specify thresholds, as well as orientation (horizontal, vertical or both).

edge_detect <- function(image, type, orientation, threshold = NULL){
  matrix_prewitt_x =  matrix(c(-1, -1, -1, 0, 0, 0, 1, 1, 1), nrow = 3, ncol = 3)
  matrix_prewitt_y = t(matrix_prewitt_x)
  matrix_sobel_x = matrix(c(-1, -2, -1, 0, 0, 0, 1, 2, 1), nrow = 3, ncol = 3)
  matrix_sobel_y = t(matrix_sobel_x)
  
  filter_x = 0
  filter_y = 0
  
  if(type == "prewitt"){
    filter_x = matrix_prewitt_x
    filter_y = matrix_prewitt_y
  }
  else if (type == "sobel"){
    filter_x = matrix_prewitt_x
    filter_y = matrix_prewitt_y
  }
  else{
    print("Error")
  }

  if(orientation == "both"){
    gradient_x <- filter2(x = image, filter = filter_x)
    gradient_y <- filter2(x = image, filter = filter_y)
    gradient_mag <- sqrt(gradient_x^2+gradient_y^2)
  }
  else if (orientation == "vertical"){
    gradient_y <- filter2(x = image, filter = filter_y)
    gradient_mag <- sqrt(gradient_y^2)
  }
  else if (orientation == "horizontal"){
    gradient_x <- filter2(x = image, filter = filter_x)
    gradient_mag <- sqrt(gradient_x^2)
  }
  else{
    print("Error")
  }
  if(! is.null(threshold)){
    gradient_mag <- threshold_values(gradient_mag,threshold)
  }
  
  return (gradient_mag)
}

Let’s test the funtion now.

Sobel operator with both orientations and no threshold.

par(mfrow = c(1,2))
sobel_image <- edge_detect(grayImage, type = "sobel", orientation = "both", threshold = NULL)
display(sobel_image,method="raster")
sobel_image <- edge_detect(myGrayImage, type = "sobel", orientation = "both", threshold = NULL)
display(sobel_image,method="raster")

Sobel operator with horizontal orientation and no threshold.

par(mfrow = c(1,2))

sobel_image <- edge_detect(grayImage, type = "sobel", orientation = "horizontal", threshold = NULL)
display(sobel_image,method="raster")
sobel_image <- edge_detect(myGrayImage, type = "sobel", orientation = "horizontal", threshold = NULL)
display(sobel_image,method="raster")

Sobel operator with vertical orientation and no threshold.

par(mfrow = c(1,2))

sobel_image <- edge_detect(grayImage, type = "sobel", orientation = "vertical", threshold = NULL)
display(sobel_image,method="raster")
sobel_image <- edge_detect(myGrayImage, type = "sobel", orientation = "vertical", threshold = NULL)
display(sobel_image,method="raster")

Sobel operator with both orientations and threshold of 0.25.

par(mfrow = c(1,2))

sobel_image <- edge_detect(grayImage, type = "sobel", orientation = "both", threshold = 0.25)
display(sobel_image,method="raster")
sobel_image <- edge_detect(myGrayImage, type = "sobel", orientation = "both", threshold = 0.25)
display(sobel_image,method="raster")

Sobel operator with both orientations and threshold of 0.5.

par(mfrow = c(1,2))

sobel_image <- edge_detect(grayImage, type = "sobel", orientation = "both", threshold = 0.5)
display(sobel_image,method="raster")
sobel_image <- edge_detect(myGrayImage, type = "sobel", orientation = "both", threshold = 0.5)
display(sobel_image,method="raster")

Sobel operator with both orientations and threshold of 0.75.

par(mfrow = c(1,2))

sobel_image <- edge_detect(grayImage, type = "sobel", orientation = "both", threshold = 0.75)
display(sobel_image,method="raster")
sobel_image <- edge_detect(myGrayImage, type = "sobel", orientation = "both", threshold = 0.75)
display(sobel_image,method="raster")

Prewitt operator with both orientations and no threshold.

par(mfrow = c(1,2))

prewitt_image <- edge_detect(grayImage, type = "prewitt", orientation = "both", threshold = NULL)
display(prewitt_image,method="raster")

prewitt_image <- edge_detect(myGrayImage, type = "prewitt", orientation = "both", threshold = NULL)
display(prewitt_image,method="raster")

Prewitt operator with vertical orientations and no threshold.

par(mfrow = c(1,2))

prewitt_image <- edge_detect(grayImage, type = "prewitt", orientation = "vertical", threshold = NULL)
display(prewitt_image,method="raster")
prewitt_image <- edge_detect(myGrayImage, type = "prewitt", orientation = "vertical", threshold = NULL)
display(prewitt_image,method="raster")

Prewitt operator with horizontal orientations and no threshold.

par(mfrow = c(1,2))

prewitt_image <- edge_detect(grayImage, type = "prewitt", orientation = "horizontal", threshold = NULL)
display(prewitt_image,method="raster")
prewitt_image <- edge_detect(myGrayImage, type = "prewitt", orientation = "horizontal", threshold = NULL)
display(prewitt_image,method="raster")

Prewitt operator with both orientations and threshold of 0.25.

par(mfrow = c(1,2))

prewitt_image <- edge_detect(grayImage, type = "prewitt", orientation = "both", threshold = 0.25)
display(prewitt_image,method="raster")

prewitt_image <- edge_detect(myGrayImage, type = "prewitt", orientation = "both", threshold = 0.25)
display(prewitt_image,method="raster")

Prewitt operator with both orientations and threshold of 0.5.

par(mfrow = c(1,2))

prewitt_image <- edge_detect(grayImage, type = "prewitt", orientation = "both", threshold = 0.5)
display(prewitt_image,method="raster")

prewitt_image <- edge_detect(myGrayImage, type = "prewitt", orientation = "both", threshold = 0.5)
display(prewitt_image,method="raster")

Prewitt operator with both orientations and threshold of 0.75.

par(mfrow = c(1,2))

prewitt_image <- edge_detect(grayImage, type = "prewitt", orientation = "both", threshold = 0.75)
display(prewitt_image,method="raster")

prewitt_image <- edge_detect(myGrayImage, type = "prewitt", orientation = "both", threshold = 0.75)
display(prewitt_image,method="raster")

We can notice both the Prewitt and Sobel operators do a very good job.
I would even say they performed better than the gradient magnitude approach because with higher threshold some values still existed.

11. Laplacian operator and Zero Crossings

We will use the Laplacian operator to calculae the second derivatives.

I will try a mask in the form: [0,1,0; 1,0,1; 0,1,0] that could be used to detect zero crossings.

zero_crossing = function(image){
  laplace_operator <- matrix(c(0, 1, 0, 1, -4, 1, 0, 1, 0), nrow = 3, ncol = 3)
  filter_operator <- matrix(c(0, 1, 0, 1, 0, 1, 0, 1, 0), nrow = 3, ncol = 3)
  
  zero_crossings_raw <- filter2(x = image, filter = laplace_operator)
  zero_crossing_raw2 <- filter2(x = zero_crossings_raw, filter = filter_operator)

  return(zero_crossing_raw2)
}

We can apply this function now to our images.

par(mfrow = c(1,2))
display(grayImage,method="raster")
display(zero_crossing(grayImage),method="raster")

par(mfrow = c(1,2))
display(myGrayImage,method="raster")
display(zero_crossing(myGrayImage),method="raster")

The results aren’t that good.

But we could improve the results using a Laplacian of Gaussian where Gaussian filtering is used to smooth out the transition and the Laplacian operator is used for edge detection.
We will use the same mask.

zero_crossing = function(image){
  laplace_operator <- matrix(c(
    0,0,-1,0,0,0,-1,-2,-1,0,-1,-2,16,-2,-1,0,-1,-2,-1,0,0,0,-1,0,0
  ), nrow = 5, ncol = 5)
  filter_operator <- matrix(c(0, 1, 0, 1, 0, 1, 0, 1, 0), nrow = 3, ncol = 3)
  
  zero_crossings_raw <- filter2(x = image, filter = laplace_operator)
  zero_crossing_raw2 <- filter2(x = zero_crossings_raw, filter = filter_operator)

  return(zero_crossing_raw2)
}

We can apply this function now to our images.

par(mfrow = c(1,2))
display(grayImage,method="raster")
display(zero_crossing(grayImage),method="raster")

par(mfrow = c(1,2))
display(myGrayImage,method="raster")
display(zero_crossing(myGrayImage),method="raster")

I think this approach with Laplacian of Gaussian gave pretty good results.